Bose-Einstein Condensation Temperature of a Homogeneous Weakly Interacting Bose 

Gas : PIMC study 
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Using a finite-temperature Path Integral Monte Carlo simulation (PIMC) method and finite- 
size scaling, we have investigated the interaction-induced shift of the phase transition temperature 
for Bose-Einstein condensation of homogeneous weakly interacting Bose gases in three dimensions, 
which is given by a proposed analytical expression T c = T < ?{1 + cian 1 / 3 + [c' 2 ln(an 1//3 ) + c^'j^n 2 / 3 + 
0(a 3 n)}, where T° is the critical temperature for an ideal gas, a is the s-wave scattering length, and 
n is the number density. We have used smaller number densities and more time slices than in the 
previous PIMC simulations [Gruter et ah, Phys. Rev. Lett. 79, 3549 (1997)] in order to understand 
the difference in the value of the coefficient ci between their results and the (apparently) other 
reliable results in the literature. Our results show that {{T c — T^)/T^}/{an^^ 3 ) depends strongly on 
the interaction strength an 1 / 3 while the previous PIMC results are considerably flatter and smaller 
than our results. We obtain ci = 1.32 ± 0.14, in agreement with results from recent Monte Carlo 
methods of three-dimensional 0(2) scalar (f> 4 field theory and variational perturbation theory. 

PACS numbers: 05.30. Jp,03. 75. Hh,02.70.Ss 



I. INTRODUCTION 

Since the recent achievement of experimental observa- 
tion of Bose-Einstein Condensation (BEC) in magneti- 
cally trapped atomic vapors yj, these inhomogeneous sys- 
tems have attracted considerable interest from experi- 
mental and theoretical sides. In dilute or weakly interact- 
ing gases, the Gross-Pitaevskii (GP) theory, a mean-field 
theory, for the condensate wave function has been enor- 
mously successful in describing the extraordinary prop- 
erties of the condensate. (For a recent review of nu- 
merous successful applications of mean field theories in 
Bose-Einstein condensation in atomic gases see Ref. (2).) 
However, the determination of the effect of repulsive in- 
teractions on the transition temperature of a homoge- 
neous dilute Bose gas at a fixed density has had a l ong 
and controversial history]! BJ1JE EJj B E El ME I 
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IMl2a|26, 27, 28, 29]. 
This non-trivial problem has been treated by different 
methods with different results. One of the reasons for 
the multitude of results and methods stems from the fact 
that the phase transition is second order, and perturba- 
tion theory typically breaks down for physical quantities 
sensitive to the collective long-wavelength modes close 
to the transition due to infrared (IR) divergences. More- 
over, the interaction, no matter how small, changes the 
universality class of the phase transition from the Gaus- 
sian complex-field model to that of the XY-mode\. 

In the weak interaction (or dilute) limit, the strength of 
the interatomic interactions can be characterized by the 
s-wave scattering length a. The shift AT C = (T e - T c °) of 
the BEC transition temperature of a homogeneous Bose 
gas away from its ideal-gas value 



behaves parametrically as 
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where m is the particle mass, ks is the Boltzmann con- 
stant, n is the number density, £(s) is the Riemann zeta 
function, £(§) w 2.612, and ci is a numerical coeffi- 
cient. Additionally, M. Holzmann et al. ]vj^ have ar- 
gued that a logarithmic term appears at order-a 2 in Eq. 
(J2J, and they have shown that this term is of the form 



h^an 1 / 3 ). They also estimated the value of the 
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numerical coefficient c! 2 using large- N s arguments in the 
three-dimensional 0(N S ) field theory. Later, Arnold et 
«L|17| were able to show that the transition temperature 
for a dilute, homogeneous, 3D Bose gas can be expressed 
to leading order as 
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where c' 2 — 19.7518 can be calculated perturbatively, 
whereas c\ and c' 2 ' require non-perturbative techniques. 
With the help of Monte Carlo data they estimated c' 2 ' « 
75.7; however, a strong debate concerns the values of the 
other numerical coefficients, especially c\ which has been 
computed using various non-perturbative methods and 
Monte Carlo techniques. There are numerous estimates 
for the parameter c\ describing the leading deviation of 
the BEC temperature due to a small repulsive interac- 
tion, and the range of variation of different predictions 
is from 0.34 (Ref. 6) to 4.66 (Ref. 21); these have been 
summarized for instance in Refs. (10, 16, 30). It ap- 
pears that the most reliable results so far are obtained 
by Monte Carlo (MC) simulations of three-dimensional 
0(2) scalar 4 field theory, c\ — 1.29 ± 0.05 determined 
by Kashurnikov et aZ.l24| and c\ = 1.32 ± 0.02 deter- 
mined by Arnold et aL[lfi|. and by variational pertur- 
bation theory (VPT), ci = 1.27 ± 0.11 determined by 
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Kastening|29j. In particular, Gruter et al.p] investigated 
the dependence of AT C numerically using Path Integral 
Monte Carlo methods. They obtained the value of c% — 
0.34 ± 0.06 after numerical extrapolation of the calcu- 
lation to the limit a — > 0. This value is about 4 times 
smaller than the above results. Holzmann et al.\lj% ar- 
gued that the difference in these results is attributable to 
a nonanalytic structure, a 2 In a, of the transition temper- 
ature in a, and that this correction gives rise to a strong 
dependence on a, even in the very dilute limit. Arnold et 
al. 0| believed that one likely problem with the PIMC 
simulation of Gruter et al.^ is inadequate system size. 
Kashurnikov et a/.|24| stated that the PIMC simulation 
did not reach the universality region because the minimal 
value of na 3 was 10 ~ 5 . 

In this paper, we attempt to understand the differ- 
ence in the value of the coefficient C\ between results 
from PIMC and the above results. We have used a 
finite-temperature path-integral Monte Carlo (PIMC) 
method 31] to study the dependence on interaction 
strength of the phase transition temperature for Bose- 
Einstein Condensation of weakly interacting Bose gases 
in three dimension. Although our PIMC method is the 
same as that used by Gruter et al., we have used smaller 
number densities and more time slices than in the previ- 
ous PIMC simulations, as motivated above. 



II. SIMULATION METHOD AND DEFINITION 
OF THE PHYSICAL QUANTITIES 

We wish to study the problem of a quantum iV-particle 
system in order to compute the phase transition temper- 
ature T c for dilute or weakly interacting Bose gases in 
three dimensions. We assume that the interparticle in- 
teraction can be described by a positive scattering length 
a, equivalent to the interaction of hard spheres of diam- 
eter a (see, e.g., Refs. [2,3]). The Hamiltonian for this 
system may be written as 

h 2 N 

* = -^ + l>. ( 4 ) 

i— 1 i<j 

where v(r) is the hard-sphere potential defined by 

v(r) = +oo (r < a) 

= (r > a). 

The statistical mechanics of quantum systems is gov- 
erned by the density matrix. For a system of N bosons at 
an inverse temperature /3, the Bose-symmetrized density 
matrix is given by 

PB (R,R'; / 3) = -i y ^p(R,PR';/3), (5) 



where R and R' are two configurations of N hard spheres. 
P denotes a permutation of particle labels among hard 
spheres, and PR is one such permutation. Evaluating 
the density matrix for interacting systems at very low 
temperatures is complicated by the fact that the kinetic 
and potential terms in the exponent of the density matrix 
cannot be separated. We can aviod this problem by in- 
serting M — 1 intermediate configurations into Eq. (J5J to 
obtain the path-integral formulation of the density ma- 
trix: 



p(R, PR; 0) = I ■■■ I dRidR 2 • • ■ dR M -i 

x p(R, Ri; r) • ■ • p(R M _ 1) PR'; r), (6) 

where f3 = l/fc^T and t = j3/M is the imaginary time 
step. The problem of evaluating the density matrix at a 
low temperature /3 _1 has been replaced by the problem of 
multiple integration of density matrices at a higher tem- 
perature t -1 . The sum over permutations in Eq.(0 com- 
bined with the integration in Eq.© can be evaluated in 
path-integral Monte Carlo (PIMC) by a stochastic sam- 
pling of the discrete paths {R, Ri,R2,- ■ ■,Raj-i,FR') 
using multilevel Monte Carlo sampling 31], an exten- 
sion of the standard Metropolis method|3^]. The PIMC 
method is essentially exact, the only necessary input be- 
ing the interparticle potential or equivalently the scatter- 
ing length a. In order to use Monte Carlo sampling, we 
must first provide a pair-product form of the exact two- 
body density matrices for the high-temperature density 
matrices that appear in the integrand Eq.@). We used 
the high-temperature approximation for the hard-sphere 
propagator derived by Cao and Berne |33|. In order to 
choose the value of M, we performed consistency checks 
by varying M to see that the results have converged. We 
used up to fifteen time-slices (M = 15). We performed 
PIMC in the continuum. The particles were confined to 
a cubic box with volume V and edge length L = V 1 ^ 3 , 
to which periodic boundary conditions were applied. We 
employed the canonical ensemble, i.e. in each simula- 
tion we fixed the temperature T, the number of particles 
N, and the dimensions L of the simulation cell. 20,000 - 
30,000 MC steps are required for equilibration depending 
on the number density. Statistical averages are collected 
from 30,000 MC steps after this. In each MC step, we 
attempted 50 - 800 trial moves at each time slice. We per- 
formed these calculations partly on IBM SP2 and partly 
on PC at the University of Georgia. The smallest-density 
simulations took about 350 CPU hours on IBM SP2. 

The scaling functions for the condensate density and 
the superfluid density are similar and imply that in the 
(dilute) interacting 3D Bose gas condensation and su- 
perfluidity occur at precisely the same temperature. Su- 
perfluid density can be calculated using the winding 
number W for simulations that have periodic boundary 
conditions |34j. Nonzero winding numbers occur when 
particles, through a series of permutations, are permuted 
with periodic images of themselves. The winding num- 
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ber is directly related to p s , the superfluid density. The 
superfluid density is given bv|34| 



p. = m (W 2 ) 
p ~ h 2 3f3N ' 

where the winding number W is defined by 



(7) 



A' 



1 dt ' 



(8) 



In order to obtain an estimate of the superfluid tran- 
sition temperature T c , we perform a finite-size scaling 
analysis. Near the critical temperature T c , the finite- 
size behavior of the superfluid density obeys the scaled 
formEl 



p s (t,L)/p = L-*/»f(tL l /»), 



(9) 



where 7r is the critical exponent of the bulk super- 
fluid fraction, p s (t)/p ~ t* , t = (T - T c )/T c , v is 
the correlation length exponent, ~ t~'\ and the 

universal function /(tL 1 '") must be analytic for fi- 
nite argument. It is assumed that 7r = v since this 
is indicated by experiment 36], renormalization-group 
calculation 37] , and the Josephson hyperscaling relation 
7r = (d — 2)v. In obtaining an estimate of the critical tem- 
perature T c , we didnot use Eq. (JSJ directly because of 
large statistical errors in our data (see Fig 2(a)). Instead 
we fit our data for different temperatures T and several 
values N of the total number of particles to a linear form 
for/ 



N x ' 3 p s (t,N)/p 



f{tN l ' 3v ) 

f(Q) + f N 1 ^(T-T c )/T c ,(10) 



where the dimensionless iV 1 / 3 has replaced the length L. 
From the fitting, we determined four fitting parameters 
/(0), / , v, andT c . 



III. SIMULATION RESULTS AND DISCUSSION 



1.3 



1.2 



1.1 



0.9 



-0.2 



0.15 



0.146 - 



0.142 



0.138 



(a) 




*N=27 
AN=64 
>N=125 
• N=216 



-0.1 



tN 



0.1 



0.2 



(b) 




: Our simulation data 
Linear fit 



0.04 



0.08 



0.12 



1/N" 



FIG. 1: Results for the noninteracting case, (a) The scaled 
superfluid fraction N 1 ^ 3 p a (t, N) /pas a function of scaled tem- 
peratures tN 1/s ". The solid line is a linear fit to Eq. JTUJ. The 
best fit parameters are the critical temperature of T c /T®= 
1.000(1) and the correlation length exponent of v — 0.90(6). 
(b) The total energy density as a function of the number 
of particles with a straight line fit to Eq. 11411 . From the 
fit we obtain the energy density in the thermodynamic limit 

(N-KX>). 



The main goal of this paper is to determine the critical 
temperature T c of a 3D homogeneous system of hard- 
sphere bosons by path-integral Monte Carlo simulations 
and finite-size scaling. We have calculated the superfluid 
fraction p s (T, N) / p for various number densities n at a 
fixed hard-sphere diameter a. 

In order to check the program carefully, we determined 
the critical temperature T® and the total energy den- 
sity in the noninteracting case, both of which are known 
exactly. Fig. 1(a) shows the scaled superfluid fraction 
N 1 / 3 p s (t, N)/ p, scaled temperarure tN 1 / 3 " ', and the re- 
sulting linear fit (solid line) to Eq. i|l(J|) for the nonin- 
teracting hard spheres. The best fit parameters are the 



critical temperature of T c /T®= 1.000(1) and the correla- 
tion length exponent of v = 0.90(6), in agreement with 
the theoretical exact values, namely T c /T® — 1 and v = 
1. 

From the expression for the specific heat c(t, L), 

dE(t,L) 



c(t, L) = 



dT 



(11) 



we obtain the energy density E(t, L), via integration, up 
to a constant 



E(t, L) = c(0, oo)T + lS a - x M"T c D{tL i -l v ), (12) 
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where t= (T — T c )/T c . Defining dD{x)/dx = g(x), we 
write the specific heat as follows : 

c(t,L) = c(0,oo)+L a ^g(tL^ u ), (13) 

where g(x) is a universal scaling function for the specific 
heat. For t — > we obtain 

E{0.L) = E c + E 1 L ( ~ a - 1) l". (14) 
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FIG. 2: Determination of the critical temperature T c at a 
fixed number density n = 5 x 10~ 3 for the interacting case, (a) 
The scaled superfluid fraction N 1 ^ 3 p s (t, N) / p as a function of 
scaled temperatures T/T°. The four lines cross at the critical 
temperature, (b) Result of our linear fit (solid straight line) to 
the data from Eq. HUH . Our estimated critical temperature 
is T c /T c ° = 1.078(1). 

The result of the fit of the energy density data to Eq. 
(|14|l is shown in Fig. 1(b). In the fit we used a = — 1 
and v — 1 for an ideal Bose gas. If we exclude the data 
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FIG. 3: Scaled superfluid fraction N 1/3 p a (t, N)/p vs tN 1/Sv 
at the smallest number density n = 1 X 10 -7 we used in 
this study. Our estimated critical temperature is T c /T° = 
1.0061(14). 



corresponding to the smallest lattice, we obtain an esti- 
mate for the energy density in the thermodynamic limit 
(N — > oo ) which agrees with the exactly known value 
within an errorbar. 

In the interacting case, we have calculated the su- 
perfluid fraction and estimated the critical temperatures 
for various hard-sphere diameters a at a fixed number 
density. We can see that the critical temperature ap- 
proaches T c = T c ° as the hard-sphere diameter decreases 
(not shown), as expected. 

In order to study the effect of interactions on the tran- 
sition temparature, we calculated the superfluid frac- 
tion p s (t, N) I p and determined the critical temperatures 
T c (n) for various number densities 1 x 1CP 7 < n < 
5 x 1(T 3 . 

Gruter, Ceperley, and Laloe@ investigated the depen- 
dence of AT C numerically using PIMC (to be referred to 
here as PRL97). In order to compare our estimated crit- 
ical temperature with PRL97 data, we have calculated 
the superfluid fraction p s /p{T,N) for N = 27, 64, 125, 
and 216 as a function of temperature T/T° at a fixed 
number density n = 5 x 10~ 3 . Fig. 2 shows our data 
(Fig. 2(a)) and our result of a linear fit (solid straight 
line) to the data (Fig. 2(b)) using Eq. lflTj|) . Our cal- 
culated superfluid fractions p s (t,N)/p are significantly 
larger than PRL97 data (cf. Fig. 1 of Ref. 6) and our es- 
timated critical temperature is T c /T® — 1.078(1), which 
is higher than PRL97 data, T c /T c ° = 1.057(2) (cf. Fig. 2 
of Ref. 6). We don't understand the reason for this large 
difference. In Fig. 3 we also show the scaled data and 
the resulting linear fit at the smallest number density n 
= 1 x 10~ 7 we used in this study. Our estimated critical 
temperature is T c /T c ° = 1.0061(14). 

The resulting dependence of the transition tempera- 
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FIG. 4: Our calculated dependence of the interaction-induced 
shift of the transition temperature of a dilute homogeneous 
Bose gas on the scattering length an 1/3 \ (a) (AT c /T°)/an 1/3 
(b) T c /T® . Whenever not shown, the error bars in our simula- 
tion data are smaller than the symbol sizes. We also plotted 
the PRL97 data (Open circles) and Eq. J3J using Arnold's 
resultsfuH (Solid line), (a) (AT c /T°)/an 1/3 depends on the 
interaction strength an}/ 3 while the PRL97 data are consid- 
erably flatter and smaller than our results, (b) The linear 
behavior is seen only at small value of an 1 / 3 . The dashed 
straight line is a fit to our four smallest data points. We 
found that the slope (ci) is 1.32(14). 



ture (a) (AT C /T®) / an 1 / 3 0) T c /T c ° on an 1 / 3 is shown in 
Fig. 4. Whenever not shown, the error bars in our sim- 
ulation data are smaller than the symbol sizes. We also 
show comparisons between our results (Filled circles) , the 
PRL97 data (Open circles), and Eq. © using Arnold's 
results (Solid line). We see that the (AT C /T C °)/W/ 3 
depends on the interaction strength an 1 / 3 while the 
PRL97 data are considerably natter and smaller than 
our results (see Fig. 4(a)). Our calculated data approach 
to Arnold's results. The linear behavior is seen only at 
small value of an 1 / 3 (see Fig. 4(b)). Our calculated data 
are slightly larger than Eq. © using Arnold's results|lq| 
(Solid line) in the linear region. The dashed straight 
line is a fit to our data points for the four smallest val- 
ues of n. We found that the slope (ci) is 1.32 ± 0.14. 
This is in agreement with the recent MC values of three- 
dimensional 0(2) scalar <p 4 field theory, ci = 1.29 ± 0.05 
(Ref. 24) and a = 1.32 ± 0.02 (Ref. 16), and the result 
by variational perturbation theory (VPT), c\ — 1.27 ± 
0.11 (Ref. 29). 

In summary, we have determined the interaction- 
induced shift of the phase transition temperature for 
Bose-Einstein condensation of homogeneous weakly in- 
teracting Bose gases in three dimensions using PIMC 
and finite-size scaling. Our results show that {(T c — 
T®)/T®} /(an 1 / 3 ) depends on the interaction strength 
an 1 / 3 while the previous PIMC resultsQ are consider- 
ably flatter and smaller than our results. We obtain c\ = 
1.32 ± 0.14, in agreement with results from recent Monte 
Carlo methods of three-dimensional 0(2) scalar <j) A field 
theory and variational perturbation theory |29j. 
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